Setup - load packages
library(ggplot2)
library(cowplot)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(colorBlindness)
library(scales)
library(bestNormalize)
library(Rmisc)
## Loading required package: lattice
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
##
## Attaching package: 'ggpubr'
## The following object is masked from 'package:plyr':
##
## mutate
## The following object is masked from 'package:cowplot':
##
## get_legend
Substrate composition
sub.tah <- read.csv('/Volumes/Macintosh HD/Users/nicolakriefall/Google Drive (nicfall@bu.edu)/Moorea_revisions/env_revised/substrate.csv',header=TRUE)
sub <- subset(sub.tah,Marine.Area!="Tahiti Faaa")
sub$Substrate <- as.factor(sub$Substrate)
#str(sub)
##total proportion
sub$prop.test <- (sub$proportion)/3
##rename sites & habitats
sub$Marine.Area <- as.factor(sub$Marine.Area)
levels(sub$Marine.Area)
## [1] "Maatea" "Tiahura"
levels(sub$Marine.Area) <- c("Mo'orea SE","Mo'orea NW")
levels(sub$Marine.Area)
## [1] "Mo'orea SE" "Mo'orea NW"
sub$Marine.Area <- factor(sub$Marine.Area, levels=c("Mo'orea NW","Mo'orea SE")) #reorder
sub$Habitat <- as.factor(sub$Habitat)
levels(sub$Habitat)
## [1] "Barrier reef" "Outer slope"
levels(sub$Habitat) <- c("Back","Fore")
levels(sub$Habitat)
## [1] "Back" "Fore"
##only interested in 2013 - year of collection
sub13 <- subset(sub,Year=="2013")
Nutrients
nuts <- read.csv("/Volumes/Macintosh HD/Users/nicolakriefall/Google Drive (nicfall@bu.edu)/Moorea_revisions/env_revised/nuts_13-14_renamed.csv",header=T)
str(nuts)
## 'data.frame': 48 obs. of 11 variables:
## $ year : int 2013 2013 2013 2013 2013 2013 2013 2013 2013 2013 ...
## $ month : int 1 2 3 4 5 6 7 8 9 10 ...
## $ date : chr "1/20/13" "2/20/13" "3/20/13" "4/20/13" ...
## $ site : chr "Tiahura" "Tiahura" "Tiahura" "Tiahura" ...
## $ habitat : chr "Fore reef" "Fore reef" "Fore reef" "Fore reef" ...
## $ replicate : int 3 3 3 3 3 3 3 3 3 3 ...
## $ Phosphate_P2O5_uM: num 0.527 0.47 0.258 0.293 0.468 0.156 0.031 0.175 0.142 0.319 ...
## $ Nitrates_NO3_uM : num 0.108 0.101 0.114 0.164 0.023 0.186 0.039 0.049 0.009 0.011 ...
## $ Nitrites_NO2_uM : num 0.069 0.008 0.028 0.036 0.013 0.006 0.019 0.002 0.034 0.005 ...
## $ Silice_SiO2_uM : num 2.58 7.11 1.28 1.08 1.2 ...
## $ Ammonium_NH4_uM : num 4.637 2.198 1.06 0.473 0.55 ...
From revisions - by month
Phosphate revised
nuts$ym <- paste0(nuts$year,"_",nuts$month)
nuts$date2 <- as.Date(nuts$date,"%m/%d/%y")
str(nuts)
## 'data.frame': 48 obs. of 13 variables:
## $ year : int 2013 2013 2013 2013 2013 2013 2013 2013 2013 2013 ...
## $ month : int 1 2 3 4 5 6 7 8 9 10 ...
## $ date : chr "1/20/13" "2/20/13" "3/20/13" "4/20/13" ...
## $ site : chr "Tiahura" "Tiahura" "Tiahura" "Tiahura" ...
## $ habitat : chr "Fore reef" "Fore reef" "Fore reef" "Fore reef" ...
## $ replicate : int 3 3 3 3 3 3 3 3 3 3 ...
## $ Phosphate_P2O5_uM: num 0.527 0.47 0.258 0.293 0.468 0.156 0.031 0.175 0.142 0.319 ...
## $ Nitrates_NO3_uM : num 0.108 0.101 0.114 0.164 0.023 0.186 0.039 0.049 0.009 0.011 ...
## $ Nitrites_NO2_uM : num 0.069 0.008 0.028 0.036 0.013 0.006 0.019 0.002 0.034 0.005 ...
## $ Silice_SiO2_uM : num 2.58 7.11 1.28 1.08 1.2 ...
## $ Ammonium_NH4_uM : num 4.637 2.198 1.06 0.473 0.55 ...
## $ ym : chr "2013_1" "2013_2" "2013_3" "2013_4" ...
## $ date2 : Date, format: "2013-01-20" "2013-02-20" ...
gg.new.pho <- ggplot(nuts,aes(x=date2,y=Phosphate_P2O5_uM,color=habitat,group=habitat,shape=habitat))+
geom_line()+
theme_cowplot()+
geom_point()+
geom_vline(xintercept=as.numeric(as.Date("2013-07-26")),linetype=4)+
geom_vline(xintercept=as.numeric(as.Date("2013-08-09")),linetype=4)+
scale_shape_manual(values=c(16,15))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
#theme(legend.position="none")+
xlab("")+
ylab("Phosphates (µM)")+
labs(color="Habitat",shape="Habitat")
gg.new.pho

Nitrates revised
gg.new.nia <- ggplot(nuts,aes(x=date2,y=Nitrates_NO3_uM,color=habitat,group=habitat,shape=habitat))+
geom_line()+
theme_cowplot()+
geom_point()+
geom_vline(xintercept=as.numeric(as.Date("2013-07-26")),linetype=4)+
geom_vline(xintercept=as.numeric(as.Date("2013-08-09")),linetype=4)+
scale_shape_manual(values=c(16,15))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
#theme(legend.position="none")+
xlab("")+
ylab("Nitrates (µM)")+
labs(color="Habitat",shape="Habitat")
gg.new.nia

Nitrites revised
gg.new.nii <- ggplot(nuts,aes(x=date2,y=Nitrites_NO2_uM,color=habitat,group=habitat,shape=habitat))+
geom_line()+
theme_cowplot()+
geom_point()+
geom_vline(xintercept=as.numeric(as.Date("2013-07-26")),linetype=4)+
geom_vline(xintercept=as.numeric(as.Date("2013-08-09")),linetype=4)+
scale_shape_manual(values=c(16,15))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
#theme(legend.position="none")+
xlab("")+
ylab("Nitrites (µM)")+
labs(color="Habitat",shape="Habitat")
gg.new.nii

Silica revised
gg.new.sil <- ggplot(nuts,aes(x=date2,y=Silice_SiO2_uM,color=habitat,group=habitat,shape=habitat))+
geom_line()+
theme_cowplot()+
geom_point()+
geom_vline(xintercept=as.numeric(as.Date("2013-07-26")),linetype=4)+
geom_vline(xintercept=as.numeric(as.Date("2013-08-09")),linetype=4)+
scale_shape_manual(values=c(16,15))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
#theme(legend.position="none")+
xlab("")+
ylab("Silica (µM)")+
labs(color="Habitat",shape="Habitat")
gg.new.sil

Ammonium revised
gg.new.amm <- ggplot(nuts,aes(x=date2,y=Ammonium_NH4_uM,color=habitat,group=habitat,shape=habitat))+
geom_line()+
theme_cowplot()+
geom_point()+
geom_vline(xintercept=as.numeric(as.Date("2013-07-26")),linetype=4)+
geom_vline(xintercept=as.numeric(as.Date("2013-08-09")),linetype=4)+
scale_shape_manual(values=c(16,15))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
#theme(legend.position="none")+
xlab("")+
ylab("Ammonium (µM)")+
labs(color="Habitat",shape="Habitat")
gg.new.amm

Results summarized - revised
ggarrange(gg.new.pho,gg.new.nia,gg.new.nii,gg.new.sil,gg.new.amm,labels="AUTO",common.legend=TRUE,legend="right",nrow=5)

#ggsave("env_nutrients_revised.pdf",height=12)
Box plots
Wrap by year
ggplot(nuts,aes(x=habitat,y=Phosphate_P2O5_uM))+
geom_boxplot()+
facet_wrap(~year)

ggplot(nuts,aes(x=habitat,y=Nitrates_NO3_uM))+
geom_boxplot()+
facet_wrap(~year)

ggplot(nuts,aes(x=habitat,y=Nitrites_NO2_uM))+
geom_boxplot()+
facet_wrap(~year)

ggplot(nuts,aes(x=habitat,y=Silice_SiO2_uM))+
geom_boxplot()+
facet_wrap(~year)

ggplot(nuts,aes(x=habitat,y=Ammonium_NH4_uM))+
geom_boxplot()+
facet_wrap(~year)

Violin plots
Phosphates
gg.pho <- ggplot(nuts,aes(x=habitat,y=Phosphate_P2O5_uM,color=habitat,shape=habitat))+
geom_violin()+
#geom_jitter(alpha=0.7)+
theme_cowplot()+
geom_boxplot(width=0.1,alpha=0)+
scale_shape_manual(values=c(16,15),labels=c("BR","FR"))+
scale_x_discrete(labels=c("BR","FR"))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
theme(legend.position="none")+
xlab("Habitat")+
ylab("Phosphates (µM)")
gg.pho

## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
nuts$habitat <- as.factor(nuts$habitat)
shapiro.test(log(nuts$Phosphate_P2O5_uM))
##
## Shapiro-Wilk normality test
##
## data: log(nuts$Phosphate_P2O5_uM)
## W = 0.95694, p-value = 0.07596
leveneTest(nuts$Phosphate_P2O5_uM,nuts$habitat)#ns
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.0922 0.7627
## 46
summary(aov(Phosphate_P2O5_uM~habitat,data=nuts)) #ns
## Df Sum Sq Mean Sq F value Pr(>F)
## habitat 1 0.0025 0.00251 0.049 0.826
## Residuals 46 2.3566 0.05123
#wilcox.test(Phosphate_P2O5_uM~habitat,data=nuts,exact=FALSE) #ns
Nitrates
gg.nia <- ggplot(nuts,aes(x=habitat,y=Nitrates_NO3_uM,color=habitat,shape=habitat))+
geom_violin()+
#geom_jitter(alpha=0.7)+
geom_boxplot(width=0.1,alpha=0)+
theme_cowplot()+
scale_shape_manual(values=c(16,15),labels=c("BR","FR"))+
scale_x_discrete(labels=c("BR","FR"))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
theme(legend.position="none")+
xlab("Habitat")+
ylab("Nitrates (µM)")
gg.nia

shapiro.test(log(nuts$Nitrates_NO3_uM)) #no
##
## Shapiro-Wilk normality test
##
## data: log(nuts$Nitrates_NO3_uM)
## W = 0.91919, p-value = 0.002784
leveneTest(nuts$Nitrates_NO3_uM,nuts$habitat)#ns
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 2.3953 0.1286
## 46
#summary(aov(Nitrates_NO3_uM~habitat,data=nuts)) #p < 0.01**
wilcox.test(Nitrates_NO3_uM~habitat,data=nuts,exact=FALSE) #p < 0.01**
##
## Wilcoxon rank sum test with continuity correction
##
## data: Nitrates_NO3_uM by habitat
## W = 435, p-value = 0.002519
## alternative hypothesis: true location shift is not equal to 0
# Wilcoxon rank sum test with continuity correction
#
# data: Nitrates_NO3_uM by habitat
# W = 435, p-value = 0.002519
sum.nia <- summarySE(data=nuts,measurevar="Nitrates_NO3_uM",groupvars="habitat")
sum.nia
## habitat N Nitrates_NO3_uM sd se ci
## 1 Back reef 24 0.171250 0.08852278 0.01806964 0.03737989
## 2 Fore reef 24 0.094625 0.05985476 0.01221780 0.02527445
Nitrites
gg.nii <- ggplot(nuts,aes(x=habitat,y=Nitrites_NO2_uM,color=habitat,shape=habitat))+
geom_violin()+
#geom_jitter(alpha=0.7)+
geom_boxplot(width=0.1,alpha=0)+
theme_cowplot()+
scale_shape_manual(values=c(16,15),labels=c("BR","FR"))+
scale_x_discrete(labels=c("BR","FR"))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
theme(legend.position="none")+
xlab("Habitat")+
ylab("Nitrites (µM)")
gg.nii

shapiro.test(nuts$Nitrites_NO2_uM)
##
## Shapiro-Wilk normality test
##
## data: nuts$Nitrites_NO2_uM
## W = 0.94972, p-value = 0.03897
#summary(aov(Nitrites_NO2_uM~habitat,data=nuts)) #p < 0.05*
wilcox.test(Nitrites_NO2_uM~habitat,data=nuts,exact=FALSE) #p < 0.01**
##
## Wilcoxon rank sum test with continuity correction
##
## data: Nitrites_NO2_uM by habitat
## W = 417, p-value = 0.008029
## alternative hypothesis: true location shift is not equal to 0
sum.nii <- summarySE(data=nuts,measurevar="Nitrites_NO2_uM",groupvars="habitat")
sum.nii
## habitat N Nitrites_NO2_uM sd se ci
## 1 Back reef 24 0.03216667 0.01396476 0.002850544 0.005896801
## 2 Fore reef 24 0.02083333 0.01636185 0.003339849 0.006909003
Silica
gg.sil <- ggplot(nuts,aes(x=habitat,y=Silice_SiO2_uM,color=habitat,shape=habitat))+
geom_violin()+
#geom_jitter(alpha=0.5)+
geom_boxplot(width=0.1,alpha=0)+
theme_cowplot()+
scale_shape_manual(values=c(16,15),labels=c("BR","FR"))+
scale_x_discrete(labels=c("BR","FR"))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
theme(legend.position="none")+
xlab("Habitat")+
ylab("Silica (µM)")
gg.sil

summary(aov(Silice_SiO2_uM~habitat,data=nuts)) #ns
## Df Sum Sq Mean Sq F value Pr(>F)
## habitat 1 0.19 0.1948 0.161 0.69
## Residuals 46 55.59 1.2084
wilcox.test(Silice_SiO2_uM~habitat,data=nuts) #ns
##
## Wilcoxon rank sum exact test
##
## data: Silice_SiO2_uM by habitat
## W = 291, p-value = 0.9593
## alternative hypothesis: true location shift is not equal to 0
Ammonium
gg.amm <- ggplot(nuts,aes(x=habitat,y=Ammonium_NH4_uM,color=habitat,shape=habitat))+
geom_violin()+
#geom_jitter(alpha=0.5)+
geom_boxplot(width=0.1,alpha=0)+
theme_cowplot()+
scale_shape_manual(values=c(16,15),labels=c("BR","FR"))+
scale_x_discrete(labels=c("BR","FR"))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"),labels=c("Back reef","Fore reef"))+
theme(legend.position="none")+
xlab("Habitat")+
ylab("Ammonium (µM)")
gg.amm

summary(aov(Ammonium_NH4_uM~habitat,data=nuts)) #ns
## Df Sum Sq Mean Sq F value Pr(>F)
## habitat 1 0.00 0.0012 0.001 0.974
## Residuals 46 50.02 1.0873
wilcox.test(Ammonium_NH4_uM~habitat,data=nuts) #ns
##
## Wilcoxon rank sum exact test
##
## data: Ammonium_NH4_uM by habitat
## W = 333, p-value = 0.3622
## alternative hypothesis: true location shift is not equal to 0
Multipanel violin plot
Results: - phosphates not sig, log-transformed anova - nitrates - nitrites - silica - ammonium
library(ggpubr)
ggarrange(gg.pho,gg.nia,gg.nii,gg.sil,gg.amm,labels="AUTO")

#ggsave("env_nutrients.pdf",width=7)
Flow
flowdat <- read.csv("/Volumes/Macintosh HD/Users/nicolakriefall/Google Drive (nicfall@bu.edu)/Moorea_revisions/env_revised/flow_rates_reformatted.csv",header=T)
ggplot(flowdat,aes(x=site,y=speed_ms,fill=site))+
geom_boxplot()+
theme_cowplot()+
scale_shape_manual(values=c(16,15))+
scale_colour_manual(values=c("#ED7953FF","#8405A7FF"))+
scale_fill_manual(values=c("#ED7953FF","#8405A7FF"))+
ylab("Flow (m/s)")+
xlab("Site-reef zone")+
theme(legend.position="none")

#ggsave("mnw.flow.pdf",width=2.5,height=2)
best.speed <- bestNormalize(flowdat$speed_ms)
flowdat$speed.t <- best.speed$x.t
a.flow <- aov(speed.t~site,data=flowdat)
summary(a.flow)
## Df Sum Sq Mean Sq F value Pr(>F)
## site 1 6128 6128 11195 <2e-16 ***
## Residuals 13536 7409 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


## hat values (leverages) are all = 0.0001477323
## and there are no factor predictors; no plot no. 5


flow.sd <- summarySE(data=flowdat,measurevar="speed_ms",groupvars="site")
flow.sd
## site N speed_ms sd se ci
## 1 MNW-B 6769 0.06481833 0.02387819 0.0002902279 0.0005689381
## 2 MNW-F 6769 0.10730170 0.02182916 0.0002653229 0.0005201163
## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] car_3.0-11 carData_3.0-4 ggpubr_0.4.0
## [4] Rmisc_1.5 plyr_1.8.6 lattice_0.20-44
## [7] bestNormalize_1.8.0 scales_1.1.1 colorBlindness_0.1.9
## [10] dplyr_1.0.7 cowplot_1.1.1 ggplot2_3.3.5
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.0 usethis_2.0.1 lubridate_1.7.10 doParallel_1.0.16
## [5] tools_4.0.2 backports_1.3.0 doRNG_1.8.2 bslib_0.2.5.1
## [9] utf8_1.2.2 R6_2.5.1 rpart_4.1-15 nortest_1.0-4
## [13] DBI_1.1.1 colorspace_2.0-2 nnet_7.3-16 withr_2.4.2
## [17] gridExtra_2.3 tidyselect_1.1.1 curl_4.3.2 compiler_4.0.2
## [21] labeling_0.4.2 sass_0.4.0 stringr_1.4.0 digest_0.6.28
## [25] foreign_0.8-81 rmarkdown_2.10 rio_0.5.27 pkgconfig_2.0.3
## [29] htmltools_0.5.2 highr_0.9 fastmap_1.1.0 rlang_0.4.12
## [33] readxl_1.3.1 farver_2.1.0 gridGraphics_0.5-1 jquerylib_0.1.4
## [37] generics_0.1.1 jsonlite_1.7.2 zip_2.2.0 magrittr_2.0.1
## [41] Matrix_1.3-4 Rcpp_1.0.7 munsell_0.5.0 fansi_0.5.0
## [45] abind_1.4-5 lifecycle_1.0.1 stringi_1.7.5 yaml_2.2.1
## [49] MASS_7.3-54 recipes_0.1.16 grid_4.0.2 parallel_4.0.2
## [53] butcher_0.1.5 forcats_0.5.1 crayon_1.4.1 haven_2.4.3
## [57] splines_4.0.2 hms_1.1.1 knitr_1.36 pillar_1.6.4
## [61] ggsignif_0.6.3 rngtools_1.5 codetools_0.2-18 glue_1.4.2
## [65] evaluate_0.14 data.table_1.14.2 vctrs_0.3.8 foreach_1.5.1
## [69] cellranger_1.1.0 gtable_0.3.0 purrr_0.3.4 tidyr_1.1.4
## [73] assertthat_0.2.1 xfun_0.27 gower_0.2.2 openxlsx_4.2.4
## [77] prodlim_2019.11.13 broom_0.7.9 rstatix_0.7.0 class_7.3-19
## [81] survival_3.2-13 timeDate_3043.102 tibble_3.1.5 iterators_1.0.13
## [85] lava_1.6.9 ellipsis_0.3.2 ipred_0.9-11